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ABSTRACT 

Context. A significant fraction of progenitors for long gamma- ray bursts (GRBs) are believed to be massive stars. The investigation of 
long GRBs therefore requires modeling the propagation of ultra-relativistic blastwaves through the circumburst medium surrounding 
massive stars. We simulate the expansion of an isotropic, adiabatic relativistic fireball into the wind-shaped medium around a massive 
GRB progenitor. The circumburst medium is composed of a realistically stratified stellar wind zone up to its termination shock, 
followed by a region of shocked wind characterized by a constant density. 

Aims. We followed the evolution of the blastwave through all its stages, including the extremely rapid acceleration up to a Lorentz 
factor 75 flow, its deceleration by interaction with stellar wind, its passage of the wind termination shock, until its propagation through 
shocked wind. 

Methods. We used the adaptive mesh refinement versatile advection code to follow the evolution of the fireball, from 3.3 seconds 
after its initial release up to more than 4.5 days beyond the burst. 

Results. We show that the acceleration from purely thermal to ultra-relativistic kinetic regimes is abrupt and produces an internally 
structured blastwave. We resolved the structure of this ultra-relativistic shell in all stages, thanks to the adaptive mesh. We comment 
on the dynamical roles played by forward and reverse shock pairs in the phase of interaction with the free stellar wind and clearly 
identify the complex shock-dominated structure created when the shell crosses the terminal shock. 

Conclusions. We show that in our model where the terminal shock is taken relatively close to the massive star, the phase of self-similar 
deceleration of Blandford-McKee type can only be produced in the constant-density, shocked wind zone. 
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1. Motivation 

There is increasing evidence that the long-duration (fcRB > 2s) 
gamma-ray burst (GRB) is associated with the collapse of a mas- 
sive star with M > 20Mo (Larsson et al. 2007). This evidence 
is supported by the association of some GRBs with a supernova 
( |Galama et al.|1998[ [Woosley & Bloo m 2006), and also by the 
association of GRBs with massive star-forming regions in distant 
galaxies (Wijers et al. 1998] [TYentham et' al. 2002). As known 
for supernova blastwave modeling, the surroundings of the ex- 
ploding stars can influence its propagation. Furthermore, some 
radio and optical observations are consistent with a scenario of 
GRB ejecta expanding into a Circ umBurst Medium (CBM) with 



wind density profile p oc r (Panaitescu & Kumar 2004 1. 



Since massive stars have significant mass-loss rates and struc- 
tured wind-blown bubbles surrounding their wind zone, we here 
investigate the ejecta dynamics as it propagates through the bub- 
ble. 

Recently, significant progress has been made in investigat- 
ing the dynamics of ultra-relativistic blastwaves expanding in 
the CBM of a Wolf-Rayet (WR) star using analytical model- 
ing ( Pe'er & Wijers |2006 1 and by numerical means, exploiting a 
Lagrangian relativistic hydro code (Nakar & Granot'2006). We 



complement these eff'orts here by a numerical simulation of the 
complete fireball dynamics, expanding in the structured CBM of 
a WR star We use grid-adaptive computations with AMRVAC 



( [Keppens et al.|2003] l to investigate the ultra-relativistic hydro- 
dynamic evolution of the fireball from its initial purely 'hot' 
phase, up to times significantly beyond its interaction with the 
transition from free-wind to shocked-wind zones. To make a 
grid-adaptive computation even feasible, this wind termination 
shock is supposed to already occur at R ~ lO'^cm, which is 
close to the progenitor compared to values given by models of 
WR evolution by [Castor et al.| ( |1973]). Since the radius of a WR 
star is within the order 3-1 IRq ( Meynet et al.|2006) , this still 
turns into a need to resolve ultra-relativistic blastwave dynamics 
over a distance of at least 6 orders of magnitude. | Van Marie et al. 



Send offprint requests to: Z. Meliani 



(2006) explored a number of physical mechanisms that could ex- 
plain a more restricted free-wind region, such as a high interstel- 
lar density and/or pressure, or a lower mass loss rate of the WR 
star One of the aims of this paper is to quantify the effect of a 
sudden, termination-shock variation in the CBM density profile 
on the dynamics of the blastwave. Some scenarios suggest re- 
producing the peculiar light-curve evolution of some GRBs (e.g. 
990123, 021211, 050904) ( [Panaitescu & Kumar|[2004l |Gendre| 
et al. 20071 by invoking an encounter of the blastwave with the 
density jump across the wind termination shock. This can lead to 
a brief brightening of the afterglow (Wijers 200 1 ). In such a sce- 
nario, the fireball expands during the first hours in a free-wind 
medium, and after several hours (or few days), the blast further 
decelerates in shock-dominated interactions with the constant 
density medium representing shocked wind ( Gendre et al. 2Q07\ . 
Before modeling the effects of such an encounter on the spectral 
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Fig. 1. The density profile of the CBM in our model, 
changes in the light curve, we here model the high energy dy- 
namics associated with this terminal shock encounter. We inves- 
tigate this for the first time numerically for blastwaves in CBM 
of massive stars. 

2. The circumburst medium 

Wolf-Rayet stars are possible progenitors GRB. During the WR 
phase, a massive star with mass M > 2OM0 has a mass-loss rate 
of Mwind - 10"^ to lO^^'Mo/yr and a fast supersonic wind with 
speeds of v„i„d = 1000 to 2500km/s ( jChiosi & MaederpW6| . 
This wind interacts strongly with the surrounding medium, cre- 
ating two shocks. At the forward shock, the ambient medium 
is swept up, where it gets compressed and heated. The reverse 
shock decelerates the wind itself, and converts almost all the ki- 
netic energy of the wind to thermal energy, producing a hot gas 



with sound speed on the order of the free-wind speed ( jWeaver 
et al. 1977| l. The two shocked regions are separated by a con- 
tact discontinuity, where Rayleigh-Taylor instabilities develop 
and lead to mixing. The shocked ISM is very dense and cools 
quickly (|Franco et al. 1994 1. However, the density in the shocked 
wind zone is lower, so its cooling time is longer. The result is an 
innermost zone with a hypersonic stellar wind and a second hot 
and almost isobaric region. It consists of shocked stellar wind, 
mixed with swept-up ISM. In fact, the total mass in this shocked 
wind region alone is dominated by mass mixed in and evapo- 
rated from the shocked ISM ahead of the contact discontinuity. 
Still, most swept-up ISM is concentrated in a third region be- 
tween contact discontinuity and forward shock, where it forms a 
thin, dense, cold shell. The mass in this thin shell remains domi- 



nant throughout the stellar-wind bubble evolution ( [Weaver et al. 
19771 1. ^ fourth region is the undisturbed interstellar medium. 

The density profile in the free stellar wind is set to nwind('") - 
Mnind/i^^v^in^r^mp). Here, is the mass of the proton. The 
shocked wind region is isobaric as the sound speed in this region 
is higher than the speed of expansion of the bubble, i.e. higher 
than the speed of the forward shock (Castor et al.|I975] Weaver 
et al.|I977 ]). In this case, the swept-up shell of ambient medium 
is driven by the pressure of the shocked wind. The pressure in 
this shocked wind region is then calculated at time t from 
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and the radius of the bubble is computed from 
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In these expressions, we introduced the density of the surround- 
ing ISM pisM- 

The transition from free to shocked wind is at the location 
of the reverse shock, and this is given by the balance between 



the thermal pressure in the shocked wind region and the wind 
ram pressure (true in the case of the energy driven phase when 
the shocked wind region cools slowly). The strong shock condi- 
tion at the reverse shock 3v^jj^^pvi'ind,eq/4 - Peq supposes that the 
pressure in the shocked wind region is much higher than in the 
far free wind. This leads to 
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In the following, we model the fireball expanding through 
only the free-wind and part of the shocked wind region. Our 
wind termination shock will be at about lO'^cm, in accor- 
dance with Eq. [3] for a high-density surrounding medium with 
PiSM = mplO^cm"-', a short WR lifetime fwR = lOOyear, a 
mass loss in the wind of M^ind = lO^^Mo/yr, and wind speed 
^wind = lO^km/s. While the last two are typical values in some 
circumburst media for GRB ( Arthur|2006) , we took a very short 
WR lifetime to investigate the case where the terminal shock is 
at a much shorter distance than the typical Ipc (Chevalier et al. 



2004). The density profile of the CBM in our model is given in 
Fig|T] The pressure profile in the free supersonic stellar wind is 
deduced by considering a Mach number M - i. The wind is 
then cold, and its pressure will have no effect on the propaga- 
tion of the fireball. In the shocked wind, the constant pressure 
is given by Eq|T] while the constant density is increased four- 
fold in accordance with the strong shock requirement. The wind 
velocity is constant throughout the wind zone and neglected in 
the shocked wind region. This is a proxy for the reduction by a 
factor of 4 expected at the termination shock. 

3. GRB shell model and evolution 

Initially we set a uniform static and hot shell that extended up to 
radius Rq - 10"cm. The initial specific enthalpy, normalized to 
c^, in the uniform shell is ?7sh = 100, and its energy is 



fish = 10^ 'ergs = 47r/?o577shC^nsh»Jp 



(4) 



where 6 stands for the thickness of the shell. We set approx- 
imately (5 ~ cAf ~ 10"cm, where A? » 3s is the dura- 
tion of the GRB. The mass of the shell is Msh - Esh/ifishC^)- 
The shell is static and hot, and the initial pressure is set to 
Psh = ^«sh (fish - 1) mp and the density from Eq.|4] Initially, 
the energy of the shell is only thermal. We use a constant poly- 
tropic index F = 4/3, as the temperature of the shell is relativistic 
and the interaction shell-ISM will be dominated by the forward 
shock, where the temperature of the shocked ISM is also rela- 
tivistic. The computation is done on a domain extending from 
a radius of lO^cm to 1.2 x lO'^cm, with the initial shell region 
between 2 x lO'^cm and 10"cm. In the grid-adaptive numerical 
simulation, we use 1200 grid points on the base grid level, and 
allow for 15 grid levels in total, with a doubling of the effective 
resolution between each grid level. 

In Fig|2] we draw the variation in time of the maximum 
Lorentz factor of the fireball. We also draw the variation of the 
Lorentz factor at the forward shock alone: this coincides mostly 
with the instantaneous maximum Lorentz factor, except for those 
intervals where reverse shock dynamics is particularly promi- 
nent, as discussed below. In Fig[3] we draw the variation of the 
maximum pressure with time. Note that we translated to comov- 
ing (at the forward shock) time for the latter. 

In the first phase of the simulated fireball, the initial hot 
shell and a fraction of swept-up matter are accelerated ther- 
mally (Fig|3]l very fast, such that a maximal Lorentz factor of 
Tmax - 75 is reached within face ~ 1200s (Fig|2]and|4|i. As the 



Z. Meliam\ R. Keppens^ '-^: GRB blastwaves through winds 



3 



shell is initially uniform, the center and back of the shell get de- 
layed in their acceleration with respect to the front of the shell, 
and they reach a lower Lorentz factor and introduce a tail struc- 
ture (Fig.|4|i. The maximum Lorentz factor reached is low in the 
sense that ymax < '7.sh = 100, due to the influence of accreted 
mass from the wind. This influence will generally depend on the 
initial energy in the shell E^h and the mass loss rate in the wind. 
In the simulation, in this acceleration phase the shell has accreted 
and accelerates with it a mass of Mace, wind - 0.024Msh. If the re- 
alized energy in the GRB is higher, the maximum Lorentz factor 
reached by the shell will increase, too. 

After this rapid acceleration phase, the swept-up wind mass 
increases enough to influence the dynamics of the shell. The 
shell now decelerates by transferring kinetic energy at both a for- 
ward and a reverse shock pair. An intermediate contact discon- 
tinuity separates shell from the swept-up, shocked wind matter. 
At the forward shock, the kinetic energy of the shell is passed to 
a swept-up shocked wind. As the shell has an internal structure 
produced during the acceleration phase (Fig. |4]i, the maximum 
Lorentz factor seen in the evolution (Fig.|2| decreases smoothly. 
This is contrasts with the sudden change seen in simulations 
where a uniform shell travels through constant medium density, 
as discussed in detail by 'Meliani et al. (20071. In this more re- 
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alistic case, the reverse shock will cross that frontal part of the 
shell with the highest Lorentz factor at about time f ~ 6 x lO^s. 
In a second phase, the overall maximal Lorentz factor remains 
constant y ~ 58 until f ~ 19 x lO-'s. This coincides with those 
times when the reverse shock crosses the middle part of the shell, 
which has acquired a constant Lorentz factor. 

After this, the reverse shock starts to propagate in the tail of 
the shell and, in this phase, the reverse shock is mildly relativis- 
tic. This lasts in Fig|2] until t ~ 9.4 hours marking a relatively 
fast decrease in the maximum Lorentz factor, which actually co- 
incides with the Lorentz factor in the, as yet, unshocked shell 
part ahead of the reverse shock. Beyond this rapid phase, the 
maximum Lorentz factor seen is situated at the front of the for- 
ward shock. In this phase more energy starts to be transferred to 
swept-up, shocked wind matter, and the Lorentz factor at the for- 
ward shock decreases with time with a slope smaller than -1/2. 
The latter dependence would be expected from by the analytic 
Blandford-McKee solution, but in our simulation not all shell 
energy has already been transferred to accreted matter, which is 
assumed in the self-similar solution. 

When the shock-dominated shell structure finally reaches the 
wind termination shock /?Rs,wind, all the matter originally in the 
free-wind zone is swept-up by the shell and compressed be- 
tween its forward shock and contact discontinuity (Fig. |4|. As 
the forward shock is relativistic (Fig. |4|, the compression ra- 
tio of the number density there is ~ 59.4 ~ (4-y(r) + 3), 
where n2 is the density in the compressed wind matter, while 
«! is the original wind density. Similarly, the thermal energy 
density in the compressed wind matter is 62 ~ (72 ~ 1)P2C^ 



(Blandford <& McKee 1976 1, which corresponds to analytical es- 



timates. However, as stated above, the shell does not yet reach 
a self-similar Blandford-McKee phase while traversing the free- 
wind zone, as the reverse shock continues to propagate in the 
tail of the shell and not all its energy is transferred to swept-up 
matter. The Blandford-McKee phase can be reached in a free 
stellar-wind region when the terminal shock would be faraway. 
Due to our assumption of a relatively close termination shock, 
our numerical simulation allows us to investigate blast waves 
interacting with wind termination shocks before the Blandford- 
McKee phase. 




Fig. 2. The variation in the maximum Lorentz factor. We also 
plot the Lorentz factor at the forward shock. The variation is 
presented a function of the distance and in the lab-frame time. 
Note the change from log to linear scale at 9 x lO'^cm. 
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Fig. 3. The variation in the (logarithm of) maximum pressure as 
a function of the distance and the comoving frame time. 

As soon as the shell structure has encountered the wind ter- 
mination shock and starts to travel through the shocked wind re- 
gion, the entire shell is actually composed of multiple regions 
(Fig [4]l. From right to left, we find (1) a forward shock now 
propagating in the shocked wind region. We then encounter (2) 
a contact discontinuity separating the re-shocked shocked wind 
from the swept-up free wind. This wind is shocked once more 
by (3) a new reverse shock propagating through the already 
swept-up compressed wind matter This matter extends to (4) 
the original contact discontinuity between shell and swept-up 
wind matter. At this new reverse shock (3), the Lorentz fac- 
tor falls to 73 ~ 10 = 0.725y2 and the density increases to 
n3 ~ 28 X lO^cm"^ = V3n2- Furthermore, the thermal energy 
63 ~ 2.1e2 and this new reverse shock remains Newtonian until 
the time t - 4.16day. The relations between the values of the 
density and Lorentz factor are in good agreement with recent 
analytical estimates by Pe'er & Wijers (2006V Between (3), the 



new reverse shock, and (4), the old contact discontinuity, there is 
at first still a region of the swept-up wind matter with values for 
n2{r) (decreased due to spherical expansion) and 72 as discussed 
above. Further leftward of (4) the old contact discontinuity, we 
find all the initial shell material, separated by (5) the old reverse 
shock, which continues to propagate in the tail of the initial shell. 
Hence, there are two regions in this part of the structure as well, 
namely shocked-shell matter between the contact discontinuity 
(4) and the initial reverse shock (5), and the original unshocked- 
shell matter 

The new reverse shock (3) eventually crosses all the region 
that extends up to the original contact discontinuity (4), and it 
arrives there at a distance Ri ~ 1.1 x lO'^cm. For times up to 
this arrival at R[, the maximum Lorentz factor seen in Fig |2] un- 
dergoes a weak variation, since the Lorentz factor 72, density 
n2, and energy in the swept-up wind region vary little. After ar- 
rival at the maximum Lorentz factor drawn in Fig |2] repre- 
sents the value in the region of shocked-shell matter ahead of 
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the new reverse (Newtonian) shock, where the Lorentz factor 
decreases more strongly. After R > 1130 x lO'-^cm, the maxi- 
mum Lorentz factor drawn in Fig|2] starts to represent the value 
at the forward shock (1) again, in the phase before it reaches the 
Blandford-McKee phase. From then on, while the new reverse 
shock continues to propagate through the original shell struc- 
ture, the maximum Lorentz factor at the forward shock starts to 
decrease with distance with a slope -0.9. This is still less than 
the slope in the expected long-term Blandford-McKee phase. To 
arrive at this phase, a still larger computational domain and a 
larger simulation time is needed. 

In Fig.|2] the time delay in the drop of the maximum Lorentz 
factor from when the shell crosses the wind terminal shock is 
thus related to the time for the new reverse shock to propagate 
through the region (2)-(4). Note that the maximal pressure in 
Fig. |3]directly increases at the wind terminal shock encounter. 

4. Conclusions 

In this work, we have investigated all phases of the GRB in a 
fireball modeled. We model the propagation of the thermal fire- 
ball through a complex wind-shaped CBM of a massive star. The 
fireball interacts with free stellar wind, crosses the wind termi- 
nal shock, and then propagates in the shocked wind zone. We 
discussed initial acceleration, energy transfer to CBM, deceler- 
ation, and interaction with the terminal shock. When the shell 
reaches the wind termination shock, the structure of the shell 
changes and multiple regions form: a forward shock, two contact 
discontinuities, and two reverse shocks characterize the evolving 
structure. We have shown that our simulations agree in terms of 
compression ratios, Lorentz factors, and energies reached at all 
these shocks with analytical estimates exploiting the jump con- 
ditions ( |Pe'er & Wijers|2006| l. 

There are pronounced differences in the deceleration as 
quantified by the variation in the instantaneous maximum 
Lorentz factor with models of the afterglow phase alone (Meliani] 
|et al.|2007) . We showed that the rapid thermal acceleration pro- 
duces a cold shell with internal structure where the Lorentz 
factor and density decrease from the head to the back of the 
shell. This has distinct consequences for its further long-term 
evolution and deceleration. In the deceleration phase before 
Blandford-McKee, the shell decelerates with a slope -0.9 in 
the constant density medium. These results obtained using high- 
resolution numerical simulation bring out important differences 
with analytical estimates, in particular with respect to the inter- 
nal structure of the expanding high-energy shell. In future work, 
the results obtained with our model will be used to deduce light 
curves, and to show how sudden rebrightenings may help to de- 
duce the position of the terminal shock. We have shown here 
that it is not appropriate to use the self-similar Blandford-McKee 
solution around the terminal shock. Future work will explore 
multi -dimensional scenarios and quantify the spectral outcome. 
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